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Summary.  In  this  project  we  have  made  strides  towards  developing  a  mesoscopic  statistical 
framework  (MSF)  to  model  non-equilibrium  ensembles  comprised  of  large  numbers  of  sim¬ 
ple  components.  As  a  testbed,  we  have  considered  inhomogeneity  formation  in  driven  non- 
equilibrium  systems,  namely  we  seek  to  predict  whether  inhomogeneities  appear  and  to  develop 
a  description  of  inhomogeneity  formation  that  is  valid  across  particle  and  system  scales.  These 
systems  can  be  tricky  to  analyze  as  mesoscale  structures  can  lead  to  highly  nonlinear  macro¬ 
scopic  responses  e.g.  near  a  critical  concentration,  the  effective  yield  stress  and  viscosity  of  a 
suspension  may  increase  by  many  orders  of  magnitude  with  the  addition  of  a  small  number  of 
particles.  Historically  these  systems  have  been  modeled  either  by  tracking  individual  particles 
(which  rapidly  becomes  prohibitively  computationally  expensive  as  the  number  of  particles  in¬ 
creases),  or  by  taking  a  phenomenological  continuum  approach  in  which  the  nonlinearities  are 
measured  experimentally  for  specific  systems  rather  than  derived  from  microscopic  principles. 

In  this  project  we  have  developed  a  hybrid  approach  which  operates  at  the  mesoscale  and 
hence  reduces  computational  cost  (relative  to  MD  approaches)  while  maintaining  a  direct  link 
between  effective  material  properties  and  microscopic  dynamics.  We  are  particularly  inter¬ 
ested  in  systems  in  which  the  microscopic  dynamics  of  the  system  are  known  and  reversible 
(e.g.  suspensions  of  particles  in  Stokes  flow).  These  assumptions  allow  us  to  take  advantage 
of  Markovian  approaches  hence  the  primary  challenges  of  the  project  centered  on  developing 
coarse-grained  statistical  descriptions  of  states  that  are  valid  at  different  length  scales,  and 
qualifying  transition  probabilities  between  states. 

Figure  1  presents  a  graphical  overview  of  our  method  which  required  the  development  of  novel 
approaches  to:  (1)  identify  clusters  and  states,  (2)  quantify  transition  probabilities,  and  (3) 
map  these  states  and  transitions  onto  a  Markovian  model.  To  focus  our  efforts,  we  chose  to 
study  the  process  of  cluster  development  in  systems  composed  of  multiple  particles  interacting 
in  simple  ways  (e.g.  dense  particulate  suspensions).  In  particular  we  have  applied  our  approach 
to  three  model  systems: 

1.  One  -dimensional  “ blocks  in  a  box. ;;  Blocks,  initially  randomly  placed  on  the  grid,  undergo 
random  motion  due  to  the  grid’s  coupling  to  a  thermal  reservoir. 

2.  Driven  spherical  particles  without  hydrodynamic  interactions.  The  particles  are  divided 
into  two  species  distinguished  by  the  direction  of  the  forcing.  Particles  do  not  experience 
hydrodynamic  forces  due  to  the  motion  of  other  particles,  but  do  undergo  thermal  Brow¬ 
nian  motion.  The  particles  interact  with  one  another  solely  through  hard-core  repulsive 
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ID 

Toy  model 

2D/3D  without 
hydrodynamics 

2D/3D  with 
hydrodynamics 

Description 

•  •• 

••• 

••• 

MSF  Transitions 

••o 

•oo 

ooo 

Markov  chain 

••• 

••• 

••• 

Why  do  clusters  form? 

••• 

•  oo 

•oo 

How  do  clusters  form? 

••o 

•  oo 

ooo 

Table  1:  Progress  to  date.  •••  =  completed,  mo  =  substantial  progress,  too  =  preliminary 
results.  MSF  =  mesoscopic  statistical  framework.  Note  that  our  focus  has  been  on  developing 
the  MSF  (white  region  of  the  table)  and  we  are  currently  finalizing  measurements  of  transition 
probabilities  in  all  three  systems.  Grey  regions  indicate  promising  new  directions  and  questions 
to  address  using  our  new  hybrid  approach. 


interactions. 

3.  Spherical  particles  with  hydrodynamic  interactions.  Particles  are  neutrally  buoyant,  sus¬ 
pended  in  a  Newtonian  fluid  and  driven  by  an  externally  imposed  shear  flow. 


In  all  three  systems  our  goal  is  to  apply  our  novel  coarse-grained  approach  to  predict  the 
conditions  under  which  clustering  occurs.  This  perspective  can  also  be  used  to  shed  light 
on  questions  related  to  why  clusters  form,  what  features  of  the  microscopic  dynamics  drive 
clustering,  and  quantification  of  growth  and  evolution  of  cluster  size  at  different  length  scales 
as  the  system  evolves.  While  we  have  made  significant  advances  along  these  lines,  our  results 
have  also  opened  up  a  number  of  new  promising  avenues  of  investigation.  Progress  of  the  project 
to  date  is  summarized  in  Table  1. 

Key  advances: 

•  Adopting  a  dendrogram  approach  to  describe  large  particle  ensembles.  This  was  an  impor¬ 
tant  insight  that  enabled  us  to  rigorously  distinguish  clustered  and  non-clustered  states. 

•  State  identification  based  on  dendrogram  properties  which  provides  a  novel  means  of 
coarse  graining  via  statistical  properties  of  the  ensemble. 

•  Quantification  of  transition  probabilities  between  states  in  three  canonical  systems  of 
increasing  physical  complexity. 

•  Integration  of  the  three  advances  above  into  a  Markovian  model  (in  progress). 
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Products: 

•  Zimoch,  McKinley,  and  Hosoi.  “Capillary  Breakup  of  Discontinuously  Rate  Thickening 
Suspensions.”  PRL  111,  036001  (2013). 

•  Zimoch  and  Hosoi.  “Cluster  identification  in  particle  suspensions.”  (In  progress). 

•  Zimoch  and  Hosoi.  “Hybrid  statistical  approaches  to  predict  clustering  in  many-body 
systems.”  (In  progress). 


Personnel:  Pawel  Zimoch  (graduate  student,  MIT),  A.  E.  Hosoi  (PI,  MIT). 
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Figure  1:  Graphical  overview  of  our  mesoscale  approach.  Coarse-grained  description,  state 
identification,  and  Markov  chain  analysis  are  all  generalizable  to  a  wide  array  of  systems.  Tran¬ 
sition  probabilities  are  unique  to  each  specific  case  and  reflect  the  physical  dynamics  relevant 
for  particular  systems. 
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A.  Background.  The  fundamental  assumption  of  equilibrium  thermodynamics  dictates  that 
all  microstates  of  an  isolated  system  can  be  observed  with  equal  probability.  In  the  majority 
of  systems  made  up  of  simple  particles,  this  assumption  leads  to  a  homogeneous  distribution 
of  particles  and  orientations.  Notable  exceptions  are  systems  in  which  topological  constraints 
make  more  phase  space  available  to  non-homogeneous  configurations,  e.g.  cubes  in  a  confined 
space,  or  elongated  rods.  In  this  project,  however,  we  focus  on  systems  made  up  of  simple 
spheres  in  which  homogeneous  behavior  is  extremely  common.  And  yet,  even  for  these  simple 
cases,  when  the  system  is  driven  away  from  equilibrium  by  an  externally  applied  force,  mass  may 
no  longer  be  homogeneously  distributed  and  clumps  or  “clusters”  may  appear.  Such  systems 
are  challenging  to  analyze  because  they  do  not  obey  the  symmetries  utilized  in  the  analysis  of 
equilibrium  systems  or  systems  that  can  be  treated  as  a  continuum. 

Historically  these  types  of  systems  have  been  modeled  using  two  types  of  strategies  (see  first  and 
second  row  in  Figure  1).  At  the  particle  scale  (top  row),  two  commonly  made  assumptions  are 
that  particles  either  do  not  interact  with  each  other,  or  that  for  each  particle  this  interaction 
can  be  averaged  over  all  other  particles,  as  if  each  particle  was  placed  in  a  “mean-field.”  In 
non-equilibrium  systems  such  as  particulate  suspensions,  where  multi-particle  interactions  can 
play  a  key  role  in  the  nonlinear  macroscopic  response  of  the  system,  the  complexity  of  particle- 
particle  interactions  quickly  renders  analysis  intractable.  In  the  current  state-of-the-art,  we 
have  a  good  understanding  of  two-particle  interactions.  Although  the  effect  of  interactions 
between  sets  of  three  particles  has  been  investigated  (see  e.g.  work  by  Jeff  Morris),  no  general 
framework  exists. 

At  the  other  end  of  the  spectrum  at  the  system  scale  (second  row  in  in  Figure  1),  equilibrium 
or  near-equilibrium  systems  can  be  approximated  as  continua,  and  long-wavelength  (hydrody¬ 
namic)  descriptions  can  be  obtained  by  averaging  on  a  scale  larger  than  individual  particles 
but  smaller  than  the  characteristic  scale  associated  with  variations  of  macroscopic  properties. 
In  non-equilibrium  systems,  such  as  polymer  melts  and  other  complex  fluids,  there  are  few 
first-principles  derivations  of  macroscopic  material  behavior.  Instead,  material  properties  are 
extracted  by  fitting  parameter  values  in  phenomenological  models  to  experimental  data. 

Bringing  these  microscopic  and  macroscopic  worlds  together  presents  a  significant  challenge. 
What  happens  at  the  intermediate  mesoscales?  In  equilibrium  statistical  mechanics,  this  ques¬ 
tion  is  answered  by  the  theory  of  renormalization,  which  shows  how  interactions  between  dif¬ 
ferent  components  of  a  system  change  as  the  descriptive  length  scale  changes.  No  equivalent 
theory  exists  in  the  non-equilibrium  world.  This  is  problematic  as  the  process  of  inhomogeneity 
development  occurs  primarily  at  the  mesoscale.  Therefore,  there  is  a  need  for  a  framework  that 
captures  the  structure  of  non-equilibrium  systems  as  they  evolve  at  scales  larger  than  a  single 
particle,  but  smaller  than  the  entire  system.  This  research  project  aims  to  address  this  gap  in 
our  understanding. 

B.  Identifying  and  coarse-graining  states.  There  are  four  essential  components  to  our  ap¬ 
proach:  1)  state  identification,  2)  coarse-graining,  3)  estimating  transition  probabilities  and  4) 
adapting  Markovian  approaches  to  capture  out-of-equilibrium  dynamics  and  to  estimate  prob- 
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Figure  2:  Illustration  of  the  dendro¬ 
gram  representation.  The  rectangle 
on  the  left  shows  the  evolution  of 
blocks  on  a  grid,  at  T  =  4,  from 
an  initial  condition  where  all  blocks 
are  stacked  at  the  center  of  the  grid. 
Color  indicates  the  number  of  blocks 
are  each  grid  point,  with  red  indi¬ 
cating  high,  and  blue  indicating  low. 
This  simulation  was  carried  out  with 
30  blocks  on  a  grid  with  100  loca¬ 
tions.  After  approximately  25, 000 
time  steps,  the  large  center  cluster 
falls  apart  due  to  the  relatively  high 
temperature.  Dendrograms  represent 
four  snapshots  at  different  times  as  in¬ 
dicated  by  the  dash  white  lines. 


abilities  of  observing  equilibrium  configurations.  The  two  historical  approaches  -  molecular 
dynamics  and  continuum  approaches  -  produce  outputs  that  are  tied  to  a  global  reference 
frame  such  as  particle  position  x(£)  or  spatially  varying  densities  p(x,  t).  In  contrast,  our  ap¬ 
proach  takes  a  state-centric  view  in  which  we  compute  the  probability  of  observing  a  state  Si 
independent  of  global  position  or  orientation  (e.g.  in  a  non  homogeneous  state,  we  typically  do 
not  need  to  know  the  positions  of  each  of  the  individual  clusters,  rather  we  are  interested  in 
the  characteristic  size  and  distribution  of  clusters). 

To  achieve  this,  we  borrow  a  strategy  from  computational  biology  that  captures  the  hierarchical 
nature  of  clustering,  namely  the  dendrogram  (see  Figure  2).  Connections  in  the  dendrogram 
represent  particles  that  are  “close”  as  defined  by  a  distance  metric  which  typically  combines 
instantaneous  relative  velocity  and  position.  States  are  defined  by  the  structure  of  the  den¬ 
drogram.  In  addition  the  to  the  appealing  state-centric  nature  of  this  type  of  representation, 
dendrograms  also  provide  an  intuitive  means  of  blending  the  discrete  nature  of  the  particles 
(typically  Lagrangian)  with  the  continuum  nature  of  the  surrounding  fluid  (typically  Eulerian). 
Finally  the  dendrograms  naturally  lend  themselves  to  coarse- graining  which  can  be  accom¬ 
plished  by  selecting  a  characteristic  scale  and  selecting  a  low  dimensional  representation  of 
detailed  structure  below  this  cut-off  (e.g.  characteristic  mean  and  variance  of  the  number  of 
particles  in  each  cluster). 

C.  Transition  probabilities:  three  model  systems.  The  dendrogram  coarse-graining  approach 
and  Markovian  dynamics  can  be  universally  applied  to  a  wide  variety  of  systems,  however  the 
physical  details  that  are  unique  to  each  individual  systems  must  be  captured  in  the  transition 
probabilities  between  states.  Ultimately  we  would  like  to  be  able  to  derive  these  transition 
probabilities  from  microscopic  dynamics  however,  to  demonstrate  the  effectiveness  of  our  ap- 
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proach,  we  begin  by  measuring  these  transition  probabilities  in  reduced  scale  simulations,  three 
of  which  are  described  below. 


Toy  model:  Blocks  in  a  box.  The  simplest  system  we  use  to  test  our  framework  is  a  one¬ 
dimensional  discrete  system  with  well-defined  transition  probabilities  between  microstates. 
Consider  a  ID  domain  with  a  periodic  array  of  N  discrete  spaces  coupled  to  a  thermal  bath 
at  temperature  T.  In  these  spaces,  there  reside  n  blocks  of  unit  mass.  There  is  no  limit  to 
the  number  of  blocks  that  can  reside  in  a  single  space;  blocks  that  occupy  the  same  space  are 
stacked  on  top  of  one  another.  There  are  no  energetic  interactions  between  the  blocks.  For 
a  block  to  pass  another  block,  the  two  must  first  form  an  intermediate  state  that  consists  of 
a  tower  of  height  h  —  2  blocks.  In  the  absence  of  block  interactions,  each  configuration  has 
the  same  energy,  so  the  blocks  undergo  a  diffusive  random  motion  on  the  periodic  array.  As  a 
result,  each  configuration  of  the  blocks  is  equally  likely.  However,  when  a  “gravitational  field” 
is  present  -  namely  an  external  potential  that  varies  with  the  vertical  position  of  the  block  - 
a  block  can  change  its  energy  via  stacking.  While  stacking  decreases  the  energy  of  the  system, 
there  are  fewer  stacked  configurations  than  unstacked  configurations,  so  entropy  drives  the  sys¬ 
tem  towards  disorder.  Thus,  the  blocks  evolve  under  two  competing  effects:  the  energy  drives 
the  system  towards  stacked  blocks,  i.e.  vertical  clustering,  and  the  system’s  entropy  drives  it 
towards  disaggregation.  The  relative  magnitude  of  these  two  effects  -  and  hence  the  equilibrium 
distribution  -  depends  on  the  temperature  of  the  thermal  bath  to  which  the  system  is  coupled. 


The  dynamics  of  the  simulation  can  take  several  different  forms,  and  the  detailed  statistics 
of  cluster  formation  depend  on  the  specific  form  of  the  chosen  transition  probabilities.  As  a 
concrete  example,  consider  Gibbs  sampling.  In  this  case,  the  system  evolves  in  discrete  time 
steps.  At  each  step,  one  particle  is  chosen  at  random,  and  a  move  to  the  left  or  to  the  right  is 
proposed  with  equal  probability.  The  energy  of  a  configuration  is  calculated  as 


E, 


' proposed 


(i) 


where  G  is  the  strength  of  the  gravitational  field  and  hi  is  the  location  above  the  base  level  of 
the  ith  block.  The  proposed  motion  is  then  accepted  with  probability 


P  = 


exp  (  Epr  oposed/T) 

exp  (Epr0p0seci  /  T^)  +  ex.p(Ecurrent/T) 


(2) 


Since  there  are  no  energetic  interactions  between  particles  themselves,  the  strength  of  the 
gravitational  field  G  can  be  set  to  1,  without  loss  of  generality.  The  behavior  of  the  system 
is  then  controlled  solely  by  the  temperature  of  the  bath,  T  which  can  be  tuned  to  shift  the 
equilibrium  distribution  from  a  clustered  to  an  unclustered  state. 

Driven  particles  without  hydrodynamic  interactions.  The  second  model  consists  of  spherical 
particles  enclosed  in  a  two-dimensional  box  with  periodic  boundaries.  Particles  interact  with 
each  other  thorough  hard-core  dissipative  interactions  which  are  critical  to  the  entropic  nature  of 
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the  force  driving  the  clustering.  There  are  no  hydrodynamic  interactions  between  particles,  but 
each  particle  individually  behaves  as  if  it  was  immersed  in  a  fluid  at  a  finite  temperature  hence 
the  particles  undergo  Brownian  motion.  After  the  system  reaches  equilibrium  at  a  prescribed 
temperature,  a  force  is  applied  to  each  particle.  The  particles  are  divided  into  two  species  driven 
with  forces  of  equal  magnitude  but  in  opposite  directions.  This  results  in  frequent  collisions 
between  particles  traveling  in  opposite  direction,  but  less  frequent  collisions  between  particles 
of  the  same  species.  As  a  result,  the  particles  separate  into  “lanes”,  each  containing  a  single 
particle  species. 


brownian 

motion 


T  =  1 .43 
Unclustered 


T  =  0.0043 
Clustered 


Figure  3:  (A)  The  forced  particle  model  consists  of  two  species  of  particles,  forced  in  opposite 
directions.  The  particles  experience  a  viscous  drag  and  stochastic  brownian  forces  due  to  the 
surrounding  fluid  bath,  but  no  hydrodynamic  interaction  forces,  (b-d)  Typical  frames  from  the 
simulation,  in  progressive  stages  of  segregation  of  particles,  (b)  Particles  just  before  the  forcing 
is  applied,  (c)  Particles  begin  to  segregation,  (d)  Long  term  steady-state  of  the  system,  with 
particles  segregated  into  two  streams. 

The  simulations  are  carried  out  using  the  HOOMD-Blue  molecular  dynamics  package,  developed 
by  the  Glotzer  lab  at  the  University  of  Michigan.  The  forcing  scheme  and  a  typical  simulation 
result  along  with  the  concomitant  dendrograms  is  shown  in  figure  3. 

Particles  in  a  shear  field,  with  hydrodynamic  interactions.  The  third  model  is  composed  of 
spherical  particles  enclosed  in  a  two-dimensional  box  with  periodic  boundary  conditions.  The 
particles  are  embedded  in  a  viscous  fluid  and  interact  with  each  other  through  low  Reynolds 
number  hydrodynamic  forces.  Therefore,  the  motion  of  any  particle  results  in  a  net  hydrody¬ 
namic  force  exerted  on  the  other  particles.  These  interactions  are  purely  dissipative  and  long 
range  (  ~  1/r).  Energy  is  input  via  a  steady  shearing  flow  field. 

The  simulations  on  this  model  are  performed  using  a  Stokesian  Dynamics  code  developed  by 
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James  Brady  at  Caltech  and  Jim  Swan  at  MIT  (formerly  Caltech).  To  discern  the  essential  na¬ 
ture  of  particle  pair  interactions,  we  begin  by  considering  an  ensemble  of  particle-pair  collisions 
in  isolation  as  illustrated  in  Figure  4.  Starting  from  random  initial  conditions,  an  ensemble 
of  particle  pairs  was  simulated  to  establish  the  long-time  particle  pair  probability  density  dis¬ 
tribution,  which  shows  the  probability  of  observing  a  particle  pair  in  a  given  configuration. 
Particles  are  far  more  likely  to  be  observed  close  to  each  other  (indicated  by  the  bright  band 
in  Figure  4,  bottom  right)  because  the  relative  motion  of  the  particles  in  these  configurations 
is  slow.  In  contrast,  the  relative  motion  of  particles  when  they  are  far  from  each  other  is  fast, 
resulting  in  low  overall  probability  of  observing  particles  in  such  configurations.  Note  that  this 
clustering  occurs  even  though  there  are  no  attractive  forces  between  the  particles;  rather  the 
time-averaged  dynamics  are  such  that  the  particles  are  far  more  likely  to  be  observed  close  to 
each  other  than  far  away. 


Figure  4:  (a)  A  pair  of  particles  in  a  periodic  bounding  box,  subjected  to  uniform  shear,  (b) 
Schematic  of  trajectories  of  the  particle  pair  in  steady  uniform  shear,  in  the  frame  of  reference 
centered  on  one  of  the  particles.  It  is  a  well-known  phenomenon  that  some  trajectories  are 
“closed”,  meaning  that  the  particles  never  separate  from  each  other.  (Bottom)  Snapshots  of  a 
pair  collision  ensemble  simulation.  (Left)  One  of  the  particles  is  located  at  the  center  of  the 
image  (not  shown),  while  the  other  particle  is  represented  by  a  single  point,  blue  or  red.  Each 
point  represents  one  simulation.  Blue  and  red  points  are  reflections  of  each  other  about  the 
center  of  the  image.  The  white  region  at  the  center  of  the  image  shows  the  excluded  volume 
due  to  the  presence  of  a  particle  at  the  center  of  the  image.  (Right)  Time-averaged  locations 
of  the  particle  pair.  Bright  colors  indicate  regions  where  the  second  particle  is  more  likely  to 
be  observed. 

D.  Markovian  model.  Having  identified  the  coarse-grained  system  states  Sn  and  transition 
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probabilities  qi^j  between  states,  it  is  a  straightforward  exercise  to  apply  a  Markovian  anal¬ 
ysis  to  evolve  the  system  though  non-equilibrium  configurations  and  compute  probabilities  of 
observing  equilibrium  states  by  invoking  principles  of  detailed  balance.  Consider  a  trajectory 
through  state  space  denoted  by: 

S  =  (S  =  m,S1,S2,...SN,S  =  n). 

Given  transition  probabilities  q^j  which  are  dictated  by  the  dynamics,  we  seek  the  evolution  of 
the  probability  of  observing  state  i  at  time  £,  i q(£),  and  the  long-term  probability  of  observing 
state  i,  7 fl.  The  first  can  be  computed  by  iteratively  applying  the  transition  probability  matrix 
to  the  initial  state.  The  second  can  be  directly  computed  as 


where 


=  log 


f  Qm-+Si  \ 


+  Vlog 

\Qsi+1^Si/ 


+  log 
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The  final  stage  of  this  project  is  the  integration  of  all  of  the  components  described  above.  In 
each  of  the  three  systems  (blocks  in  a  box,  particles  without  hydrodynamics,  and  particles  with 
hydrodynamics)  we  can  represent  states  of  the  system  by  dendrograms  (examples  shown  in 
Figures  2  and  3).  These  dendrograms  are  then  coarse-grained  and  classified.  Finally  transition 
probabilities  between  these  classifications  can  be  measured  numerically  (currently  underway) 
and  predictions  from  the  Markov  model  can  be  compared  with  full  MD  simulations. 
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